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(DGLAP) equations for the evolution of fragmentation functions using the 
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PROGRAM SUMMARY 



Title of program: evolve 

Computer and operating system: Program tested on Sun running SunOS 5.7, Alpha running 
OSF1 4.0, Dell-PC running Linux Mandrake 7.2. 

Programming language used: C++ with g++ compiler 
No. Lines in distributed program: 2500 

Keywords: UHECR, fragmentation functions, DGLAP equations, Laguerre method, super- 
symmetry. 

Nature of Physical Problem: In order to predict the spectra of UHECR produced by the decay 
of a dark matter superheavy particle with mass Mx one needs to calculate fragmentation functions 
at the energy scale Mx- These can be calculated from measured low energy fragmentation functions 
using the DGLAP equations. 

Method of solution: The DGLAP equations are solved by expanding them in Laguerre polyno- 
mials which reduces their integration to the computation of a set of coefficients. These coefficients 
are given by algebraic recursive relations. 

Typical running time: A few seconds. 

Restriction of the program: Gluon coherence at x < 10~ 3 is not included. 
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LONG WRITE-UP 



I. INTRODUCTION 

Over one hundred cosmic ray events with an energy higher than 4 x 10 19 eV have been observed 
by different observatories (see [|J for a recent review). Were their sources at cosmological distances, 
d > 50 — 100 Mpc, they would interact with the Cosmic Microwave Background (CMB) on their 
way to the Earth and lose their enormous energy. Therefore, one expects the sources for these 
Ultra High Energy Cosmic Rays (UHECR) to be not far from the Galaxy. However, there a few 
astrophysical sites in the galactic neighbourhood which could accelerate a charged particle to such 
high energy. Faced with this conundrum it has been suggested that UHECR are not accelerated 
at all but are created at this ultra high energy by the decay of massive dark matter particles X 
generated in the early universe In order to test this hypothesis one needs to calculate the 

angular distribution of events || and the spectra ||-|f| produced by the decay of the population 
of X particles clustered in the galactic halo.. In the present work we will concentrate on the 
calculation of the spectra. 

The X decay contribution to UHECR is proportional to the inclusive decay width of particle 
X, with mass Mx, into particle h PJ 



1 dr(X -►/» + ...) /-i dz 1 dr a (y, n 2 , M- 



Tx dx „ Jx z T a dy 



x, 



D h a (z,^). (1) 

y=x/z 



Here x is the fraction of the energy of X carried by h and z is the fraction of the energy of parton a 
carried by h. The first factor in the integrand, the decay width of X into parton a, dr a /dy, is 
calculable in perturbation theory. It encapsulates all kinematical effects in many-body decay ||. 
In lowest order and for two-body decay it is proportional to 5(1 — y). The second factor, the 
non-perturbative D%, is the fragmentation function (FF) for particles of type h from partons of 
type a. The mass scale \x is the factorization scale, \i ~ Mx [flCfl - Particle h is any final state: n, 
p, 7, e, u e , Ufj, or v T . 

The FF satisfy a set of coupled integro-differential equations, the Dokshitzer-Gribov-Lipatov- 
Altarelli-Parisi (DGLAP) equations Jll| , p^ ]. Given experimental data at some low energy scale, 
say Mz, an initial set of FF D^(x, M z ) can be extracted and evolved using the DGLAP equations, 
to obtain the FF at some higher scale D^(x, M x ). 

Hadronic structure functions (SF) satisfy as well DGLAP equations. Although DGLAP equa- 
tions for SF are similar to DGLAP equations for FF they are not equivalent. Several approaches 
have been taken to solve the DGLAP equations, usually applied to SF. One is the Mellin transform 



method [13 which transforms the integro-differential equations into ordinary differential equations. 
However, at the end one needs to invert the Mellin transform to find the solution in terms of x 
which is a process with notorious numerical problems. The QCDNUM program defines a 
grid in x and the energy scale /j?. The calculation of the SF on the grid points is based on the 
computation of convolution integrals that are evaluated as weighted sums. 

A very elegant method to solve the DGLAP equations was introduced by Furmahski and 
Petronzio JTf|. It expands these integro-differential equations using Laguerre polynomial so that 
their integration is reduced to calculating a set of numerical coefficients using simple algebraic 
recursive relations. The Laguerre method has been implemented in numerical codes by several 
groups |16-|lg[|. 
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Usually numerical codes deal with SF, and all of them are oriented to collider physics. Several 
include polarization and other aspects which add complexity to the numerical evolution and that 
are irrelevant for cosmic ray studies. Motivated as well by present or past experiments, most of the 
numerical codes evolve SF using Standard Model equations, without considering supersymmetry 
(SUSY) or any model beyond the Standard Model (SM). The theoretical basis of SUSY DGLAP 
evolution for SF was studied in Refs. [ p!9|j20[| . Numerical solutions for the evolution of SF in several 



supersymmetric scenarios have recently been presented in Ref. [21]. In UHECR studies, including 
SUSY is of paramount importance. Most of the evolution from the low energy scale Mz to the high 
energy scale Mx will be governed by SUSY equations as long as SUSY is a symmetry of nature 
and the SUSY breaking scale is of the order of the weak scale. 

We have written a numerical code to evolve FF using the DGLAP equations. It has been 
written bearing in mind its application to UHECR and therefore does not include physics aspects 
such as spin dependent functions that are relevant for collider physics but not for UHECR physics. 
It includes SM evolution and SUSY evolution. We have chosen the Laguerre method to solve 
numerically the DGLAP equations. We have generalised the Laguerre method to include SUSY 
evolution, which is not a trivial task, at least in the way that this method is presented in the 
literature. 

Matrix algebra becomes very important when solving SUSY DGLAP evolution by the Laguerre 
method. Therefore, it is convenient to code the algorithm using an object oriented language. We 
have chosen C++, which in the present context provides a good framework to perform matrix 
calculations. The use of templates available in this programming language simplifies the code. The 
final result is a fast algorithm with good accuracy for the relevant range of x and // 2 . 



II. DGLAP EQUATIONS 
A. General Equations, Leading Order and Notation 

The DGLAP equations can be written as 

dD Q a ^f } = E « S (^ 2 )) ® Dfrx, (2) 

n fx k 7I " 

where a s (fj, 2 ) is the strong coupling constant and Pba( x > a s) is the splitting function for the parton 
branching a — ► b. Here the convolution of two functions A(x) and B(x) is defined as 

A(x)® B(x) = f 1 — A{z)B{-). (3) 

Jx Z Z 

The splitting functions can be expanded perturbatively 

P ba (x, a s ) = P ba {x) + 0{a s ) (4) 

We will limit our study to leading order in a s and therefore ignore 0(a s ) corrections to the splitting 
functions. 

It is also convenient to define the following dimensionless evolution parameter 
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b being the coefficient in the leading order beta function governing the running of the strong 
coupling: f3{a s ) = —ba 2 s . We take D% to represent the sum of particle h and, if different, its 
antiparticle h. 

B. Standard Model Equations 

The Standard Model DGLAP equations for the evolution of fragmentations functions are well- 



known 1 10 , 22 1 . There are two parton species: quarks q^, k = 1, . . . np an d gluons g, with np the 
total number of flavours. Conventionally, one defines the following linear combinations (for ease of 
notation the superscript h is omitted) 

D qt =D qk +D qk (6) 

k 

D q7 = D qk - D qk (8) 



D Q k = D gt ~ ^D q . (9) 



The non-singlet functions D q - and Z?Q fc obey the equations 

d T D q -=P qq ®D q - (10) 
drD Qk =P qq ®D Qk , (11) 

while the evolution of the singlet function D q is coupled to that of the gluon function D g as 

*te)-(':*£-)-ft)- 



The splitting functions were calculated in Refs. pi. 22]. Given the FF at some initial scale no for 
the quarks and gluon g, Eqs. (|i~C|-12) completely determine their evolved values at some other 
scale fj,, to leading order in a s . 

C. Supersymmetric Equations 

In a supersymmetric model (SUSY), besides the quarks and gluon one has the superpartners: 
squarks Sk and gluinos A. In addition to the linear combinations (§) one now defines 

D st =D Sk +D Sk (13) 

D ^E D st ( 14 ) 

k 

D-=D Sk -D Sk (15) 

Ds k ^ D s+ - ±D S . (16) 

k Hp 

The non-singlet function D - and D - evolve together as do DQ k and D$ k 

k 
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(17) 
(18) 



The singlet functions for quarks and squarks, D q and D s , are coupled to the gluon and gluino 
functions, D g and D\, as 



(19) 
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In leading order the SUSY Eqs. (^-19) allow us to calculate the FF for all quark and squark 
flavours, gluons and gluinos at some scale [i once their values at some initial scale no are known. 
The SUSY DGLAP equations have been given in the literature to leading order for structure 
functions |T^ , pO| , Here we have presented their form for FF. It is easy to see that one just needs 
to transpose the matrix elements keeping the n F factors in the same place to move from SF to FF 
equations. The SUSY splitting functions were calculated in Refs. 



III. LAGUERRE ALGORITHM 

A. Evolution Operator 

In numerical studies it is better to consider the quantity xD(x,t). The evolution equations 
for xD are the same as those for D if one multiplies the splitting functions by x. This improves 
numerical stability since the 1/x singularity shown by splitting functions with a gluon in the final 
state cancels off. In general one has 

d T xD(x,T) = xP(x) ®xD(x,t), (20) 

where D is a (i-dimension vector and P is a d x d matrix whose elements are splitting functions. In 
the SM d = 1 for Eqs. (ID and <^ while d = 2 for Eq. @. In a SUSY model d = 2 for Eqs. @ 
and fllSl) whereas d = 4 for Eq. (p|). 

Following Furmariski and Petronzio one introduces the evolution operator E(x,t) 

xD{x, t) = E(x, t) ® xD{x, 0). (21) 

The evolution operator is a d x d matrix which satisfies the following integro-differential equation 

E{x,t) = xP{x)®E{x,t) (22) 

and the initial condition 

E(x,0) =5(l-x). (23) 
One introduces the Laguerre expansions 
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e- x P(e- x ) = J2(xP) n L n (x) (24) 



n=0 

oo 



E(e- x ) = Y,E n L n (x), (25) 
n=0 

with L n (x) the Laguerre polynomials of order n, which form an orthonormal basis in the interval 
(0, oo) with weight exp —x. Equivalently, L n (— Inx) are an orthonormal base in the interval (0, 1) 
with unity weight. A key property of the Laguerre polynomials is their closure under convolution 

/ dz L n (z)L m (x - z) = L n+m (x) - L n+m+1 (x). (26) 
Jo 



Substituting Eqs. (|24]) and ( [25|) into Eq. using Eq. (|2q ) and the fact that the Laguerre 

polynomials form a vectorial base one gets the following system of linear equations 

n 

En(r) = J2 (xP)n-mE m (T), (27) 
m=0 

being 

xP = xP , (28) 
xP m = xP m - xP m -i m > 1. (29) 



Laguerre expansion of the initial condition Eq. (|23| ) translates into 

E n {r = 0) = /. (30) 

The Laguerre expansion transform an integro-differential equation into a set of ordinary differential 
equations. Thus there is no need to perform any intricate quadrature. One is left with an infinite 
number of equations n = 0, 1 . . . oo. In practice one truncates the Laguerre expansion at some 
finite n = NLAST. 

B. Expansion of the Splitting Functions 

The coefficients of the Laguerre expansion of a function can be quickly calculated if its Mellin 
transform is known. In the present work we are interested in the product 

F{x) = xP(x), (31) 

where P(x) can be any splitting function. The Laguerre series is 

oo 

F(e-y) = J2E n L n (y). (32) 

n=0 

The Laguerre coefficients F n are given by 

/■oo 

F n = / dye-yL n (y)F(e-y), (33) 
Jo 

where L n (y) is the Laguerre polynomial of order n. Let F(s) be the Mellin transform of F{x) 
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F(s) = / dxx s ~ L F(x) 
Jo 

then one can calculate the F n 's by means of the formula [p 



1 
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oo 

£ 

n=0 



F n U T 



(34) 



(35) 



For the SM the Mellin transform of the splitting functions to leading order have been listed in many 
works, see for example [^]; for a SUSY model the Mellin transforms of the splitting functions 
to leading order are given in Ref. |jl9| . In both cases the Mellin transforms of P(x) are linear 
combinations of only six simple functions: 1, l/(s + n) with n = —1, 0, 1, 2 and a{s) = 7 + ip{s + 1) 
where 7 is the Euler constant and ip(s) is the digamma function. Using Eq. ( |35|) it is easy to show 
that each of these six simple functions that appear in the Mellin transform (xP)(s) = P(s + 1) 
gives a fixed contribution to F n . The rules to calculate the Laguerre coefficients of xP(x) from the 
Mellin transform of P(x) are then: 



As usual one defines 



P(s) 



(xP(x)) n 
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(36) 
(37) 

(38) 

(39) 

(40) 

(41) 



Z{n) = X; (-l) m d.) C(m + 1) 



(42) 



m=l 



where C( n ) is t^ 16 Riemann £ function. All the Laguerre coefficients for SM and SUSY splitting 
functions are listed in Appendix [A|. 



C. Algorithm in d dimensions 

Let us find the recursive relations that will allow us to solve the DGLAP equations for any 
arbitrary dimension d. The quantities E n and xP n are d x d matrices. In analogy with the d = 1, 2 



formulae presented in 15], let us try the following ansatz for the solution to Eq. (^ 



d n k 

^n(r)=E eAiT ElF^>> (43) 
i=l k=0 
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where the Aj's are the eigenvalues of xPq. For every i, k, n with < k < n B\ n is a constant d x d 
matrix; we have in total ^(NLAST + 1)(NLAST + 2) x d x d 2 matrix elements. Substituting 
Eq. ( f43|) into Eq. ( |27|) we obtain the following two matrix relations for the Bf n , d > 1, 

= (xP - XJ) B% n (44) 
B^ 1 = (xPo - XJ) B\ n + J2 *Pn- m Bt m , (45) 

m=k 

while substitution of Eq. (H) into Eq. (^) gives 

^ = E^,n- (46) 
i=l 

The matrix equation Eq. ( (44| ) is equivalent to (iVLj4ST + l) x <i x d(<i— 1) algebraic equations, since 
det (xP - XJ) = 0. Equation (|5|) is equivalent to ±N LAST(N LAST + 1) xdxd 2 equations and 
Eq. ( |46| ) is equivalent to (NLAST + 1) x d 2 equations. Adding up these three number we obtain 
^(NLAST + 1)(N LAST + 2) x d x d 2 algebraic relations, which is the same as the total number 
of matrix elements. Therefore, the set of equations (f44|- 46) completely determine all the matrix 
elements of Bf n . They will be given as recursive relations. In order to find them one defines the 
projectors associated with the matrix xPq 

n ^ n ^M (4?) 

■ / ■ Aj An 

They have the following properties: 

d 

E n * = / ( 48 ) 

i=l 

UiUj = UiSij (49) 
(xP - Xil) ILj = (Xj - Ai)IL; (50) 

As the key step we define suitable auxiliary matrices for d > 1: 

bin = (51) 
bl n = Bl n -Y,{^-M) k ^Bln, 0<k<n. (52) 

Using Eqs. (^-46) and the projector properties Eqs. (f48|-|5"o|) we have found that all the b^ n and 
B\ n satisfy the following relations: 

B% = Ik (53) 

= (54) 
n-l 

= (xP Q - XJ) b\ n + ]T xPn-mB^ (55) 

m=k 

n-l 

BiX 1 = (xPo ~ A,/) B\ n + ]T *Pn- m Bl m . (57) 

m=k 
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In the particular case d = 2 these relations reduce to the ones given in The above recursive 
relations can easily be implemented in a computer using an object oriented language. Then we 
can numerically calculate the B\ n once the matrices xP n are given. The evolution operator is 
calculated using Eq. (|4^) without need to integrate numerically its differential equation (p?|) using 
any finite step algorithm. 

Finally, the evolved FF is given by the truncated series 

NLAST n 

xD(x,T)= ^ En-m(T)(xD) m 
n=0 m=0 

x (L n (-lnx) - L n+1 (-lnx)) , (58) 
where (xD) m are the coefficients in the Laguerre expansion of the initial FF 

oo 

xD(x,t = 0) = Y,( xD )™L m (-lnx). (59) 

m=0 



IV. NUMERICAL ANALYSIS 
A. Evolution Steps 

For clarity we will assume flavour universality in the decay of X, hence we will only consider 
the coupled singlet quark and gluon evolution, Eq. (O) for the Standard Model (SM), and coupled 



singlet quark, gluon, singlet squark and gluino evolution, Eq. (19) for a SUSY model. Particular 
models for X may have different branching ratios for different flavours. The source code includes 
all the routines necessary to evolve each quark and squark flavour. 

A (s)parton is not included in the evolution as long as the energy scale is lower than its mass; 
when its threshold production scale is crossed, it is added to the evolution equations with an 
initially vanishing FF and it is assumed to be a relativistic particle. 

In the SM case the code evolves the q and g initial FF from Mz to Mt, the top quark mass, 
with the number of flavours set to rap = 5, and then evolve from Mt to Mx with rip = 6. Taking 
n-p = 6 in the whole range from Mz to Mx does not introduce any significant difference in the 
final spectrum. 

In the SUSY case the code evolves the q and g initial fragmentation functions from Mz to the 
supersymmetry breaking scale Mstjsy > Mt using the SM equations to obtain D^(x, Mg UgY ), with 
i = q,g. Then it takes D^{x, M| ugY ), i = q,g, and Dj(x, M| ugY ) = 0, j = s, A, and evolves them 
from Mstjsy to Mx using the SUSY equations. All spartons are taken to be degenerate with a 
common mass -Msusy • In the context of structure functions a broken SUSY scenario with different 



masses for s and A was studied |21], showing no significant difference with models with a unique 
SUSY mass. One expects the same result to hold for FF. 

In UHECR one is interested in the final spectra of baryons p + n, neutrinos (sum of all three 
families) and photons. The evolution of neutrinos and photons is equivalent, from a numerical 
point of view, to that of baryons. In the present work we will concentrate on baryons. Initial FF 
are extracted from LEP data at the energy scale Mz (see || for further details on initial FF and 
evolution of baryons, neutrinos and photons). 
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B. Results 



Let us first show SM evolution. In Fig. [T] we plot the fragmentation functions for baryons from 
quarks and gluons at the scale Mz (fit to LEP data) and their evolved shape at Mx = 10 10 , 10 12 
and 10 14 GeV. Following the standard convention in UHECR studies we always plot the quantity 
x 3 D a (x, [i 2 ). One can see that as the final scale increases the number of particles grows at low x 
and diminishes at high x, a well-known result from many previous studies of scaling violations. 

Next we compare SM evolution of FF with SUSY evolution. In Fig. g we show the common 
initial baryon curves at Mz, their shape after SM evolution up to Mx = 10 12 GeV, and their 
evolved shape at the same final scale after SUSY has been switched on at Msusy = 400 GeV. It is 
clear that the SUSY curves have evolved further than the SM curves. The difference between the 
two scenarios stems chiefly from the different running of a s (fi 2 ). In a SUSY model a s decreases 
with the growth of energy scale more slowly than in the SM because of the increased contribution 
to the /3-function from the SUSY partners. Since the rate of change di nfl 2D a (x, ^ 2 ) is proportional 
to a s (see Eq. (|2|)), a larger a s translates into a larger amount of evolution. In other words, for 
the same initial and final scales we obtain tsusy(Mx) > tsm(Mx), using Eq. (||). 

Figure || shows the quark and gluon functions at Mz, their evolved values at Mx = 10 12 GeV 
using the SUSY equations for scales larger than Msusy and the radiatively generated squark and 
gluino functions, all at the same final scale Mx- We find that starting from vanishing values at 
-^SUSY the squark and gluino functions start to grow and catch up with the quark and gluon 
functions, respectively, at small x and from scales a few orders of magnitude higher than Msusy- 
This behaviour can be understood qualitatively if one bears in mind that at low x the leading 
splitting function for quarks is 2n-pP gq ~ AnpCp/x, which is equal to the leading splitting function 
for squarks 2npP gs ~ An-pC-p/x. For gluons and gluinos the leading splitting functions tend as well 
to a common value, P gg ~ 2Ca/x and P g \ ~ 2Ca/x, which is however different from that of quarks 
and squarks. 

SUSY evolution does not depend strongly on the chosen supersymmetry breaking scale Msusy- 
In Fig. | we show the curves obtained taking Msusy = 200, 400 GeV and 1 TeV. The higher the 
value of Msusy > the less evolved the final curves for q and g. This follows from our comparison of 
the SM evolution and the SUSY evolution. If SUSY switches on later (higher Mstjsy) the energy 
range over which the SM equations hold is larger. As we have seen already, DGLAP evolution is 
slower when just the SM equations are employed. 

C. Convergence, Accuracy and Scope 

The final FF generated by the Laguerre algorithm are accurate for large values of the evolution 
parameter r. For very small values of r the evolution operator approaches a delta function 5(1 — x) 
and thus truncation of the Laguerre series generates loss of accuracy. The final result is stable for 
intermediate values of NLAST, the values of n at which the Laguerre expansion is truncated. In 
Fig. | we show SUSY evolution of q and g to the final scale M x = 10 12 GeV, for NLAST = 9, 12, 15. 
We see that the curves approach to a common curve in most of the range of x. For large values of 
NLAST roundoff errors start to accumulate. In double precision, for SM evolution NLAST has 
to be smaller than 30 while for SUSY evolution, where the number of matrix operations increases 
substantially, NLAST has to be smaller than 20. 
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In the limit x — ► the Laguerre polynomials with n > L n (— Inx) go to infinity. One needs 
a large number of polynomials to achieve good precision. However, taking NLAST too large, 
roundoff errors become important. The code give good accuracy for x larger than a few time 10 -3 . 
In any case, for x < 10~ 3 the DGLAP equations @ no longer hold 22]. For very small x gluon 
emission coherence has to be taken into account; this modifies the kernel of the DGLAP equations. 
We have not included gluon coherence in the present version of the code. 

The functions xD(x,r) fall rather fast when x — > 1. This fact slows down the convergence 
of the Laguerre series in the large x region [jl5[| . If the asymptotic behaviour goes as xD(x) ~ 
(1 — x) a when x — * 1 then one can improve the convergence by extracting the term (1 — x) a 
explicitly. The Laguerre series is then rewritten in terms of generalised Laguerre polynomials 
L { n\-lnx) @. Alternatively, one can take an analytical approach and calculate the asymptotic 
exponent of xD(x,t) in the limit x — > 1 P]. Nonetheless, one has to keep in mind that at large x 
FF are not well measured at low energy, hence, even in the case of perfect numerical accuracy in 
the evolution, the final result would not be reliable for x > 0.6. 



V. DESCRIPTION OF THE INPUT AND OUTPUT 

To speed up the calculation the code takes as input the Laguerre coefficients of the initial xD 
at the scale Mz- These are given in the default files bar.ldat for baryons, gam.ldat for photons 
and nu.ldat for neutrinos. They have been obtained from LEP data (baryons) or from the QCD 
Montecarlo generator HERWIG (2*| (photons and neutrinos), see Ref. ||] for further details. These 
files contain two fields: the Laguerre coefficients for the quark singlet function and the Laguerre 
coefficients for the gluon function. The routine laguerre . cc is provided in case one wishes to start 



with a different set of initial data and needs to calculate the Laguerre coefficients, see Sec. VL 

The main variables of the code are set with the help of a simple user interface. There are two 
setup levels. The level setup only allows the user who runs the code to set the final scale in 
the evolution Mx which is stored in the variable double Efinal. In this case the evolution is 
supersymmetric with SUSY breaking mass Mstjsy = 400 GeV and baryons are evolved. Level 1 
gives more freedom to the user. Besides the choice of Mx, now the user has to decide between SM 
evolution or SUSY evolution (variable bool SUSY), in the case of SUSY evolution Mstjsy must be 
set (variable double SUSYMass). The user has to select one particle to evolve: baryons, photons 
or neutrinos. The variable string ParticleData stores the name of the file with initial data 
for the selected particle. The standard template library class string is included with the header 
<string>. 

If xD(x) ~ (1 — x) a , a > 0, when x — > 1, the variable alpha can be set to an integer larger 



than to improve convergence for large x, as explained in Sec. IV C. Level defaults alpha=0. 

There are other variables that are initialised to default values, which are the recommended 
values. To change them one needs to edit the main program, recompile and link. The variable 
double Einit stores the initial scale in the evolution. The default value is Mz but it could be set 
to any value with M^ < Einit < Mt. Remember to recalculate the initial Laguerre coefficients if 
Einit is modified. The variable int NLAST stores the number of term in the Laguerre expansion 
(n = 0, 1 . . . , NLAST). The variables xmin and xmax store the minimum and maximum values of 
x, respectively, for the output of xD(x) while NSTEPS controls the number of points in between. 

The code evolves xD up to the final scale Mx and writes the final functions to standard output. 
Output has three fields for the SM (x, xD q , xD g ) and five in a SUSY model (x, xD q , xD g , xD s , 
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VI. DESCRIPTION OF THE CODE 



The following classes are defined in the code: 

• User, UserO, Userl. TheyaredefinedinthefilesUser.handUser.ee. They encapsulate the 
user interface that sets up the physical parameters of the run. The polymorphic class User 
contains the virtual function setParameters () , which prints a short help to standard error. 
The classes UserO and Userl inherit class User and redefine setParameters () . The func- 
tion UserO: : setParameters () manages the level setup while Userl: : setParameters () 
manages level 1. 

• StandardModel, SUSYModel. They are defined in the files className.h and className.ee. 

They encapsulate the physical model. Their member functions calculate the evolution 
parameter r and the splitting functions P\, a . The class SUSYModel inherits the class 
StandardModel since at low energy any SUSY model must include the Standard Model. 
The class StandardModel inherits class AuxFunc (see next). 

• AuxFunc, Laguerre. They are defined in the files DtherClasses .h and DtherClasses.ee. 
The class AuxFunc defines auxiliary functions required to calculate the Laguerre expansion 
of the splitting functions. In particular it stores tabulated values of the Riemann £ function. 
The class Laguerre encapsulates the Laguerre polynomials and the generalised Laguerre 
polynomials ]l5|. 



Vector, Matrix. They are defined in Arrays. h and Arrays.cc. They are template classes 
with one template parameter which is the vector or square matrix size. The vector and 
matrix operators needed for the Laguerre algorithm are overloaded. 

xDlXl, xD2X2, xD4X4. They are defined in the header files className.h. These classes 
perform the numerical integration by calculating the matrices B using the recursive relations 



given in Sec. Ill C, The matrices B and the evolution operator E are members of these classes. 



Storage for the B matrices is optimised taking into account that some indices are triangular, 



see Eq. @). The class xDlXl integrates the SM Eqs. (|lg) and ([□]), xD2X2 integrates the 
SM Eq. ([ID and the SUSY Eqs. (|17|) and (^), and xD4X4 integrates the SUSY Eq. (0). All 
three classes are template classes with one generic data type which can be StandardModel or 
SUSYModel. All three classes inherit classes AuxFunc and Laguerre. In the present version 
of the code only the classes xD2X2<StandardModel> and xD4X4<SUSYModel> are used in 
the main program; the class xDlXl is provided for completeness. 

The main program Evolve . cc takes the Laguerre coefficients of the initial FF for q and g at 
Mz and calculates the final FF at Mx- It encodes the following steps: 

1. It prompts the user for input and initialises variables. It declares the input file stream qgFile 
(class if stream included with <f stream>) for the input quark singlet and gluon FF. 

2. It declares the objects sm of type StandardModel and susy of type SUSYModel. 
The public member functions StandardModel: : setNFlavoursAndTau(EInit ,EFinal) and 
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SUSYModel : : setTau(susyMass ,EFinal) will be called to set the number of flavours (5 below 
Mt, 6 above) and calculate the evolution parameter r. 

3. It declares as well the objects qg of type xD2X2<StandardModel> and qgsl of 
type xD4X4<SUSYModel>. The former solves the coupled evolution of singlet quark 
and gluon in the SM, while the latter solves the coupled evolution of singlet 
quark, gluon, singlet squark and gluino in a SUSY model. The public mem- 
ber functions xD2X2<StandardModel> : : setB(sm) , xD2X2<StandardModel> : : setE(sm) 
calculate the matrices B and the evolution operator E in the SM whereas 
xD4X4<SUSYModel>: :setB(susy), xD4X4<SUSYModel> : : setE(susy) calculate B and E 
in a SUSY model. 

4. The member function setlnitialxDLaguerre () in class xD2X2<StandardModel> is over- 
loaded; it accepts an ifstream type argument or void argument. It reads the initial La- 
guerre expansion coefficients of xD(x,t) or matches Laguerre coefficients at the Mt scale. 
The mutator xD4X4<SUSYModel> : : setlnitialxDLaguerre (qg) matches final SM Laguerre 
coefficients to initial SUSY Laguerre coefficients at Msusy- 

5. The public member functions xD2X2<StandardModel> : : setFinalxDLaguerreO and 
xD4X4<SUSYModel> : : setFinalxDLaguerre () calculate the evolved Laguerre coefficients. 
The generalised Laguerre coefficients are calculated with the member functions 
setFinxDGenLagO . Finally the accessor getxD(xl,xr,NSTEPS, alpha) calculates the 
evolved fragmentation functions xD(x,t) and writes the result to standard output. If 
alpha=0 the Laguerre coefficients are used, if alpha is an integer larger than the gen- 
eralised Laguerre coefficients are used instead. 

The header file CommonDefs.h sets the precision to double, defines some numerical factors 
and sets the maximum number of terms in the Laguerre expansion NMAX=30. The header file 
PhysicalConst .h defines some physical constants used in the computation. 

As mentioned in Sec. [V| default files with Laguerre coefficients for the initial FF are provided. 
If one wishes to use another set of initial data, one will need to calculate their Laguerre coefficients. 
The file laguerre . cc provides a routine to calculate these coefficients when one has xD in tabulated 
form, i.e. one has a file with two fields: the firsts one is x and the second one xD{x). The routine 
in laguerre . cc declares an object of class Laguerre and therefore must be compiled with the file 
DtherClasses.ee (the file Makefile is provided with the source files to facilitate compilation). 
Notice that the routine in laguerre.ee is independent of the main program in Evolve.cc. The 
Laguerre coefficients in Eq. (|5^) are given by 

/•oo rl 

(xD) n = dye- y (xD)(e- y )L n (y)= dxxD(x)L n (-lnx). (60) 
Jo Jo 

The integral is performed numerically using the extended trapezoidal rule. Previous implemen- 
tations of the Laguerre method tend to use parametric fits to the data of the sort Ax a (l — x) 13 
at the initial scale. In this case the Laguerre coefficients can be calculated analytically and are 
given in terms of infinite sums. However, the x range of validity of this fits tends to be rather 
narrow. Furthermore, initial data is naturally obtained in tabulated form from experiments and 
from Montecarlo generators. Therefore we prefer calculating the initial Laguerre coefficients as 
explained above without use of any intermediate parametric fit in x space. 
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VII. CONCLUSIONS 

In order to test models where UHECR are produced by the decay of superheavy dark matter 
of mass Mx one needs to calculate the predicted spectra for baryons, photons and neutrinos. The 
X produced injection spectrum for any primary cosmic ray is basically given by fragmentation 
functions at the energy scale Mx • We calculate FF by evolving low energy FF up to the ultra high 
energy Mx using the DGLAP equations. We have solved numerically the DGLAP equations using 
the Laguerre method. In our study we have considered two scenarios: Standard Model evolution 
and SUSY evolution. We have generalised the Laguerre method to include supersymmetry. The 
final result is a numerical code which is fast and accurate for 2 x 10~ 3 < x < 0.6 and Mx > 1000 
GeV. 
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APPENDIX A: LAGUERRE COEFFICIENTS FOR THE SPLITTING 

FUNCTIONS 



The Laguerre expansions of the splitting functions can be calculated using Eq. (|35|). Let us 



first list the Laguerre coefficients of xP a b(x) in the SM [15,17|: 
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Next, we present the Laguerre coefficients of xP a b(x) in a SUSY model. To calculate them we 
use the table of Mellin transform given in 1£] (in this reference there is a mismatch between the 
normalization of the splitting functions and the normalization of the DGLAP equations that we 



correct). From the rules given in Subsec. QIB we obtain 
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TEST RUN INPUT AND OUTPUT 

% evolve 

Evolve fragmentations functions up to the energy scale M_X 
Two setup levels : 
Level -> Set M_X 

Level 1 -> Set M_X, Standard Model or SUSY evolution, 

M_SUSY (in the latter case) , particle (baryon, 
photon or neutrino) , and alpha 

Enter level (0 or 1) : 1 

Enter final scale M_X in GeV: lel2 

Enter for Standard Model evolution, 1 for SUSY evolution: 1 
Enter SUSY breaking scale M_SUSY in GeV: 400 
Enter particle (b baryons, g photons, n neutrinos): n 
Enter alpha, alpha >= (type if you don't know): 4 

Fragmentations functions evolved up to M_X=le+12 GeV 
SUSY evolution with M_SUSY=400 GeV 

Initial Laguerre coefficients read from file: ./init_data/nu.ldat 
Parameter alpha=4 

0.001 413.595 40.1029 360.98 37.8363 
0.00114574 372.335 35.3804 319.938 33.179 
0.00131271 334.828 31.14 282.839 29.0114 
0.00150402 300.755 27.3395 249.372 25.291 
0.00172321 269.822 23.9402 219.247 21.9782 
0.00197435 241.756 20.9062 192.194 19.0364 
0.00226209 216.308 18.2042 167.958 16.4314 
0.00259176 193.247 15.8036 146.305 14.1317 
0.00296947 172.365 13.6762 127.012 12.1079 
0.00340223 153.468 11.7957 109.873 10.3329 
0.00389806 136.381 10.1384 94.6951 8.78167 
0.00446615 120.945 8.68196 81.2995 7.4311 
0.00511703 107.013 7.40617 69.5181 6.25986 
0.00586277 94.454 6.29235 59.1949 5.24839 
0.00671719 83.1462 5.32339 50.1848 4.37872 
0.00769614 72.9806 4.48362 42.3533 3.63442 
0.00881775 63.8575 3.75871 35.5756 3.0005 
0.0101028 55.686 3.13559 29.7366 2.46334 
0.0115752 48.3833 2.60235 24.7301 2.0106 
0.0132621 41.8736 2.14816 20.4588 1.63113 
0.0151949 36.0872 1.76324 16.8335 1.31494 
0.0174093 30.9602 1.43872 13.773 1.05307 
0.0199465 26.4336 1.16663 11.2037 0.837566 
0.0228534 22.4526 0.939831 9.05914 0.661396 
0.026184 18.9664 0.75193 7.27963 0.518366 
0.03 15.9278 0.597254 5.812 0.40307 
0.0343721 13.2927 0.470784 4.60909 0.310818 
0.0393814 11.0197 0.368103 3.62942 0.237572 
0.0451207 9.07067 0.285354 2.83672 0.17988 
0.0516964 7.40964 0.219183 2.19954 0.134819 
0.0592305 6.00339 0.166701 1.69086 0.0999292 
0.0678626 4.82113 0.125434 1.28761 0.0731663 
0.0777527 3.83453 0.0932815 0.970304 0.052841 
0.0890841 3.01766 0.0684755 0.722594 0.0375722 
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0.102067 2.34699 0.0495397 0.530891 0.0262409 
0.116942 1.80132 0.0352533 0.383979 0.0179478 
0.133985 1.36174 0.0246149 0.272665 0.0119761 
0.153511 1.01153 0.0168107 0.189457 0.00775903 
0.175883 0.736047 0.0111848 0.128271 0.00485098 
0.201516 0.522559 0.00721274 0.0841796 0.00290404 
0.230884 0.360077 0.00447843 0.0531969 0.0016483 
0.264532 0.239156 0.00265422 0.0321031 0.000876103 
0.303084 0.151691 0.00148427 0.0183037 0.000429735 
0.347255 0.0907179 0.000770882 0.00972273 0.000191815 
0.397863 0.0502393 0.000363635 0.00472289 7.78663e-05 
0.455846 0.0250863 0.00015074 0.00204541 3.02731e-05 
0.52228 0.0108333 5.21486e-05 0.000762316 1.30163e-05 
0.598395 0.00376894 1 . 37897e-05 0.000232293 6.639e-06 
0.685603 0.000921878 2.35586e-06 5.33011e-05 3.14451e-06 
0.785521 0.000115102 1.76694e-07 7.61964e-06 8.88797e-07 
0.9 2.33062e-06 1.74687e-09 2.87264e-07 5 . 27823e-08 
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FIG. 1. Standard Model fragmentation functions for baryons from quarks and gluons, at the 
initial scale M z (solid lines) and the final scales: M x = 10 10 GeV (dashed line), M x = 10 12 GeV 
(dotted line) and Mx = 10 14 GeV (dot-dashed line), showing scaling violations. 
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FIG. 2. Fragmentation functions for baryons from quarks and gluons, at the initial scale Mz 
(solid lines) and evolved to a decaying particle mass scale of 10 12 GeV, for SM evolution (dotted 
lines), and, the more pronounced, SUSY evolution (dashed lines) taking Mstjsy = 400 GeV. 
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FIG. 3. Dependence on Msusy- the dashed lines are final (s)parton functions with 
Msusy = 200 GeV, the solid lines have Msusy = 400 GeV and the dotted lines have 
-Msusy = 1 TeV. The initial and final scale are always Mz and Mx = 10 12 GeV, respectively. 
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FIG. 4. SUSY evolution with M x = 10 12 GeV and M SU sy = 400 GeV. The solid lines are 
obtained with NLAST = 15, the dotted lines with NLAST = 12 and the dashed lines with 
NLAST = 9. All curves have been calculated with alpha=0. While good convergence is achieved 
very fast at intermediate values of x, for x > 0.2 one needs to take a larger number of terms in the 
Laguerre expansion in order to damp out unphysical oscillations. 
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